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COMPUTATIONS  OF  CHAOTIC  FLOWS  IN  MICROMIXERS 


1.  Introduction 

Components  of  microdevices  containing  fluids  are  now  integral  parts  of  micro-electro-mechanical 
systems.  These  microfluidic  components  are  used  to  transport  particles  and  reactants,  and  they 
are  also  used  as  chambers  in  which  to  perform  various  types  of  chemical  and  biological  analysis.  A 
general  and  critical  problem  in  the  design  of  microfluidic  components  is  how  to  create  a  physical 
environment  in  which  the  movement  and  location  of  fluids  (gases  and  liquids)  and  entrained  par¬ 
ticulates  can  be  controlled  for  mixing,  separation,  and  delivery.  Sometimes  it  is  just  necessary  to 
move  two  separate  streams  to  some  location,  other  times  it  is  necessary  to  have  two  streams  mixed 
to  a  certain  degree  at  a  specified  site,  and  sometimes  it  is  necessary  to  deliver  a  large  particle  in 
the  flow  to  a  specified  surface  for  detection  or  modification.  General  and  useful  references  include 
Ho  (1998),  Gad-el-Hak  (1998,  2002),  Karniadakis  et  a 1.  (2002),  Nguyen  &  Wereley  (2002),  and 
Breuer  (2005). 

In  current  terminology,  a  component  is  active  when,  for  example,  there  are  additional  forces 
or  driving  mechanisms,  such  as  mechanically  driven  parts,  unsteady  electric  fields,  or  imposed 
oscillations  in  the  inflow  (see,  for  example,  Evans  et  a 1.  (1997),  Rife  et  a 1.  (2000),  Glasgow  & 
Aubry  (2003),  Tsai  &  Lin  (2002),  Ukita  &  Kanehira  (2002),  and  Yang  et  a 1.  (1998)).  Often, 
but  not  always,  this  implied  an  attempt  at  active  control  of  some  property  of  the  flow  system.  A 
passive  component  usually  is  one  whose  geometry  is  set  before  it  functions,  its  parts  do  not  move, 
and  there  are  no  additional  fields  or  forces  or  time-dependent  effects  imposed  on  the  system  (see, 
for  example,  Bessoth  et  a 1.  (1999),  Park  et  a 1.  (2004),  Hong  et  a 1.  (2004),  Johnson  et  a 1.  (2002), 
and  Stroock  et  al.  (2002)).  In  general,  the  first  attempt  is  to  design  a  passive  system  because  it  is 
simpler  to  design  and  requires  less  energy  to  operate. 

When  there  is  substantial  internal  structure  in  the  channel  or  the  component  is  passive,  the 
flow  can  become  chaotic  (see,  e.g.,  Jones  et  al.  (1989);  Khakhar  et  al.  (1987);  Ottino  (1989); 
Wiggins  &  Ottino  (2004);  Chien  et  al.  (1986);  Kaplan  et  al.  (2004,  2005)).  Here  the  adjective 
chaotic  is  not  meant  in  the  colloquial  sense  of  a  state  of  extreme  confusion  and  disorder.  Rather,  it 
is  meant  as  it  is  used  in  nonlinear  dynamics  theory  to  mean  a  dynamical  system  that  is  extremely 
sensitive  to  its  initial  conditions,  a  system  that  produces  complicated  patterns  that  are  not  really 
random  (e.g.,  Ottino  (1990);  Ottino  &  Wiggins  (2004)).  (We  will  go  more  into  this  later  in  the 
paper.)  In  this  case,  the  velocity  field  itself  can  still  be  relatively  straightforward  to  compute,  but 
mixing  and  scalar  transport  become  complex  as  a  series  of  finer  and  finer  striations  form  in  time 
M  anuscript  approved  December  8,  2005. 


or  distance.  Sometimes  it  is  possible  to  use  simplified  models  of  nonlinear  dynamics  to  analyze 
these  flows,  and  these  give  important  insights  into  the  distinct  regimes  of  the  flow  (e.g.,  Stroock  & 
McGraw  (2004);  Khakhar  et  a 1.  (1986)).  A  problem  with  totally  relying  on  these  models  is  that 
they  are  not  usually  quantitatively  correct  and  might  not  be  qualitatively  correct  in  the  presence  of 
other  physical  effects  such  as  molecular  diffusion,  thermal  conduction,  chemical  reactions,  or  large 
particles  in  the  flow.  In  this  case,  we  can  only  use  numerical  simulation  of  the  governing  equations. 

Liquids  in  microcomponents  are  characterized  by  low  flow  velocities,  constant  density,  and 
hence  low  Reynolds  numbers  (Re  «  1  —  10).  The  governing  equations  are  the  multispecies  in¬ 
compressible  Navier-Stokes  equations,  generally  including  terms  representing  molecular  diffusion, 
chemical  reactions,  and  interphase  drag  and  energy  transfer.  Thus  the  flow  is  laminar  and  often  in 
the  creeping  flow  regime.  These  equations  are  relatively  straightforward  to  solve  numerically  when 
the  chamber  geometry  is  simple  enough,  the  walls  are  relatively  smooth  or  have  well  described 
functional  behavior,  and  the  component  is  passive.  (That  is,  there  are  no  external  forces  besides 
those  which  cause  the  material  to  flow  at  a  constant  rate  through  the  chamber.)  When  these 
simplifications  are  not  valid,  obtaining  the  solution  becomes  extremely  costly  and  time  consuming. 

This  paper  describes  numerical  simulations  of  a  three-dimensional,  geometrically  complex, 
incompressible,  “creeping”  flow  that  is  a  prototype  for  one  used  in  a  micromixer.  The  original 
motivation  for  this  work  was  to  develop  an  affordable  computational  approach  that  would  solve 
the  Navier-Stokes  equations  accurately  enough  on  a  single  workstation  to  design  microfluidic  com¬ 
ponents.  Then,  combined  with  new  methods  for  rapid  prototyping  and  testing,  we  could  use  the 
numerical  model  to  determine,  a  priori,  the  best  configurations.  A  major  hurdle  for  computing 
many  of  these  very  low-Reynolds-number  flows  in  complex  geometries  is  resolving  the  effects  of 
nonlinear  dynamical  processes,  and  determining  how  these  complicate  the  flow  structure.  The 
difficulty  of  the  problem  is  equivalent  to  that  of  turbulence  in  higher-Reynolds-number  flows.  As 
with  turbulence,  there  are  a  number  of  low-order  approaches,  here  based  on  nonlinear  dynamics, 
for  describing  the  effects  of  chaotic  advection  on  mixing.  The  next-level  problem,  the  one  we  are 
facing,  is  to  understand  the  flow  mechanistically  and  to  develop  methods  to  actually  compute  it 
and  predict  the  behavior  quantitatively  as  well  as  qualitatively. 

This  paper,  a  greatly  extended  expanded  version  of  Kaplan  et  a L  (2004),  then  describes  our 
efforts  to  compute  an  incompressible  flow  that  becomes  more  and  more  complex  and  convoluted, 
developing  ever- narrowing  striations  of  interpenetrating  fluid  as  it  moves  downstream.  Producing 
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our  best  result  involved: 


-  Developing  new  numerical  algorithms  and  approaches  to  computing  incompressible  flow  that 
allow  us  to  take  advantage  of  the  accuracy  of  nonlinear  monotone  methods  for  convection  (Liu 
et  al.  2005), 

-  Developing  an  approach  to  complex  geometry  that  would  not  require  paying  the  cost  of  un¬ 
structured  grids,  and 

-  Analyzing  the  flow,  as  the  calculations  evolved,  so  that  we  could  take  advantage  of  established 
physical  properties  of  the  flow,  use  a  number  of  computational  and  theoretically  motivated 
shortcuts,  and  still  arrive  at  a  reasonable  final  solution  given  current  resources. 

Here  we  use  the  Stroock  et  a 1.  (2002)  staggered-herringbone  micromixer  as  a  guide  for  testing 
the  approach.  They  described  their  component  as  a  microchannel  containing  a  series  of  bas-relief 
herringbone  structures  (ridges  in  the  channel)  where  one  sequence  consists  of  a  total  of  twelve 
herringbone  structures,  six  aligned  one  way,  followed  by  six  aligned  oppositely.  The  channel  in  the 
original  experiment  consisted  of  15  sequences,  two  of  which  are  shown  schematically  in  Fig.  1.  This 
figure  also  indicates  the  dimensions  used  in  the  computations  described  later  in  this  paper.  Half 
of  the  liquid  entering  the  channel  contained  a  glycerol-water  mixture  (indicated  by  A )  and  half 
contained  this  same  mixture  with  a  small  amount  of  light-sensitive  fluoroscein  dye  (indicated  by 
B ).  Images  taken  by  Stroock  et  al.  (2002)  with  a  confocal  micrograph,  Fig.  2,  show  how  striations 
of  A  and  B  form  as  the  flow  moves  down  the  channel.  By  the  time  the  fluid  has  passed  the  last 
sequence,  it  is  effectively  intermixed.  Note  that  if  the  diffusion  coefficient  between  A  and  B  had 
been  higher,  better  mixing  would  have  occurred  sooner. 


2.  Description  of  the  Method 


The  physical  model  of  the  micromixer  requires  solution  of  the  incompressible,  multispecies  Navier- 
Stokes  equations, 

V-v  =  0  (1) 

^  +  V-(pvv)  =  -VP  +  /xV2v  (2) 

dc ■ 

+  (v  •  Vcj)  =  Dij  \72Ci,  i  =  l,n8  (3) 
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where  v  is  the  fluid  velocity,  p  is  the  fluid  density,  P  is  the  pressure,  p  is  the  dynamic  viscosity 
coefficient,  q  is  the  concentration  of  component  i  of  which  there  are  ns,  and  Dij  is  the  binary 
diffusion  coefficient  between  two  components  i  and  j.  These  equations  are  solved  on  a  three- 
dimensional  Cartesian  mesh  using  a  finite-volume  formulation  with  a  colocated  arrangement  of 
variables  (Liu  et  a 1.  2005).  The  advective  fluxes  are  calculated  with  flux-corrected  transport  (FCT) 
(See  Oran  &  Boris  (2001)  for  details),  using  the  package  LCPFCT  (Boris  et  a 1.  1993),  which  is  a 
a  high-order,  nonlinear,  monotone  algorithm  for  solving  continuity  equations  with  source  terms. 

2.1  Tests  of  the  Computations :  Unperturbed  MicroChannel  Flow 

First  consider  a  three-dimensional  rectangular  channel  with  the  dimensions  4000  x  90  x  200  /xm3, 
no  obstacles  in  the  flow,  and  smooth  walls.  (This  is  similar  to  that  shown  in  Fig.  1,  but  without  the 
ridges.)  The  channel  contains  water,  characterized  by  a  dynamic  viscosity  coefficient  of  0.014  g/cm- 
s.  The  inflow  velocity  is  specified  at  1  cm/s,  the  outflow  boundary  condition  is  specified  at  1  atm, 
and  the  walls  are  no-slip.  Two  passive  scalars,  labeled  A  and  B,  are  introduced  at  the  inflow. 
Scalar  A  represents  water,  B  represents  water  with  fluoroscein,  and  the  binary  diffusivity  between 
A  and  B  is  taken  as  1  x  10-6  cm2/s.  Scalars  A  and  B  follow  the  streamlines  of  the  flow  but  do 
not  interact  with  the  fluid  or  each  other.  Here  Re  ^  1. 

Figure  3  shows  the  computed  velocity  profile  in  the  x-  and  ^-directions,  pressure  distribution, 
and  concentration  of  B.  The  computational  cells  are  5  pm  in  all  directions,  so  that  there  are 
800  x  18  x  40  computational  cells  in  the  grid.  For  this  flow,  boundary  layers  form  along  the 
microchannel  walls  where  they  grow,  and  merge,  and  eventually  produce  a  parabolic  pipe  flow 
velocity  profile.  Along  the  centerline,  the  axial  velocity  accelerates  to  approximately  2  cm/s.  After 
the  start-up  region,  the  steady  state  pressure  gradient  along  the  channel  is  linear  and  is  2950 
dynes/cm3.  The  only  signficant  velocity  in  the  ^-direction  is  at  the  inflow  boundary,  due  to  the 
interaction  of  the  fluid  with  the  microchannel  walls.  At  the  inflow,  the  two  fluids  are  completely 
unmixed.  As  the  fluids  move  down  the  microchannel,  there  is  a  small  amount  of  mixing  at  the 
interface  due  to  mass  diffusion.  (This  is  seen  by  the  slight  spreading  of  the  green  area  between  the 
red  and  purple  regions.) 

Figure  4  shows  these  concentrations  at  the  inflow  and  exit  planes,  x  =  0  and  4000  /xm, 
respectively.  At  the  exit  plane,  the  components  are  better  mixed  at  the  fluid  interface  at  the  upper 
and  lower  walls,  y  =  0  and  y  =  90  pm.  This  is  because  the  fluid  velocity  is  slowed  considerably 
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along  the  walls,  so  that  the  two  fluids  have  a  longer  residence  time  there  during  which  they  can 
diffuse  and  mix. 

The  results  shown  in  Figs.  3  and  4  were  obtained  using  5  /im  x  5  /im  x  5  /im  computational 
cells.  Additional  calculations  at  grid  resolutions  of  2.5  /xm  x  2.5  /xm  x  2.5  /xm  and  10  /xm  x  10  /xm 
x  10  /xm  gave  steady-state  pressure  gradients  of  2925  and  3025  dynes/cm3,  respectively.  From  this 
essentially  quadratic  convergence,  the  steady-state  pressure  gradient  is  estimated  as  approximately 
2912  dynes/cm3.  An  analytic  solution  for  the  incompressible  Navier-Stokes  equations  for  flow  in  a 
rectangular-shaped  duct  is  given  by  White  (1991).  Using  the  conditions  appropriate  for  our  sim¬ 
ulation,  the  analytic  solution  predicts  a  steady-state  pressure  gradient  of  2895  dynes/cm3.  Figure 
5  shows  a  comparison  of  the  x-velocity  for  point  values  of  the  analytic  solution  and  cellaveraged 
values  in  the  simulation  for  a  grid  resolution  of  5  /xm.  The  profile  shows  x-velocity  along  the  y  and 
z  axes,  in  the  region  of  fully  developed  flow. 

2.2  Tests  of  the  Computations:  Effects  of  Staircased  Grid 

Because  we  are  using  a  Cartesian  grid  and  the  herringbones  are  set  at  angles  to  the  grid,  we  need 
to  quantify  the  effects  of  staircasing  on  the  solution.  As  an  example  of  the  tests  we  have  performed, 
consider  the  problem  shown  in  Fig.  6,  flow  over  a  triangular-shaped  obstacle.  The  computational 
domain  consists  of  100  x  20  x  40  cells,  where  each  cell  is  5  /xm,  corresponding  to  a  domain  of  500  /xm 
x  100  /xm  x  200  /xm.  The  obstacle  is  located  in  the  region  200  <  x  <  250  /xm  and  y  <  35  /xm, 
and  the  staircasing  is  in  the  z-direction.  In  one  test  case,  the  obstacle  staircasing  is  at  the  grid 
resolution  (5  /xm),  and  in  another  case,  the  staircasing  is  coarser,  at  10  /xm  resolution. 

Figure  7  shows  results  for  the  two  cases  at  the  y  =  25  /xm  plane.  Because  this  is  such  a 
highly  viscous  flow,  there  is  essentially  no  difference  in  the  two  test  cases.  Figure  8  shows  line 
profiles  from  this  plane  at  z  —  100  /xm  (along  the  centerline).  The  flat  region  in  each  profile, 
where  215  <  x  <  250  /xm,  corresponds  to  the  obstacle  location.  For  these  two  cases,  there  are  no 
differences  in  x-velocity,  ^/-velocity  and  pressure,  and  only  minor  differences  in  z-velocity.  From 
these  and  other  similar  tests  (not  shown),  we  conclude  that  with  enough  resolution,  we  should  be 
able  to  use  staircasing  for  these  computations. 


5 


3.  Computation  of  the  Herringbone  Mixer  Flow 

Channel  flow  is  characterized  by  a  pressure  drop  from  the  entrance  to  the  exit.  This  is  a  function  of 
channel  length,  cross-section,  viscosity  and  flow  rate.  The  total  pressure  drop  along  the  length  of 
the  channel  increases  when  there  are  obstacles  present.  The  pressure  drop  is  linear  along  the  length 
of  a  channel  for  an  incompressible  fluid.  Therefore,  if  we  do  not  have  any  obstacles  and  we  are 
given  the  pressure  drop  (dp/dx)  across  some  segment  of  the  channel,  we  then  know  the  pressure 
drop  across  any  segment.  Below  we  show  and  use  the  fact  that  given  a  repeating  sequence  of 
obstacles,  as  we  have  with  the  herringbone  patterns  in  Fig.  1,  we  can  make  an  analogous  statement 
about  selected  segments  of  the  channel.  Then  we  use  this  information  to  help  decrease  the  cost  of 
a  many-segment  calculation. 

For  idealized  scalar  transport,  where  there  is  no  feedback  from  the  scalars  and  the  background 
flow,  we  can  separate  the  problem  into  two  types  of  computations.  First,  we  compute  the  steady 
velocity  field  by  solving  the  Navier-Stokes  equations,  and  then  we  use  this  to  compute  the  motion 
of  the  scalar  field  subjected  to  this  velocity.  By  first  obtaining  the  velocity  field,  we  can  vary  how 
(that  is,  location,  time  history)  passive  scalars  are  injected.  Knowing  the  velocity  field  is  the  key 
to  understanding  the  flow  in  a  geometry. 

3.1  The  Velocity  Field 

The  inflow  conditions  and  the  dimensions  of  the  herringbone  mixer,  as  shown  in  Fig.  1,  are  the 
same  as  those  of  the  microchannel  with  no  obstacles.  The  individual  herringbones  are  20  /xm  in 
height,  spaced  150  /xm  apart,  and  their  surfaces  in  the  flow  are  no-slip  walls.  The  channel  contains 
two  sequences  of  herringbone  structures,  where  each  sequence  consists  of  six  herringbones  with  the 
short  segment  on  the  left  and  the  long  segment  on  the  right,  followed  by  six  more  herringbones 
with  the  short  segment  on  the  right  and  the  long  segment  on  the  left.  The  ratio  of  the  length 
of  the  long  segment  to  short  segment  is  2:1.  Results  of  grid  resolution  tests  will  be  discussed 
throughout  this  paper,  however,  the  velocity-field  computations  discussed  in  detail  below  are  from 
an  intermediate-resolution  computation  in  which  Ax  =  Ay  =  A z  —  5  /xm,  requiring  800  x  18  x  40 
computational  cells  for  two  sequences. 

The  major  effect  of  the  herringbone  structures  is  to  deflect  the  flow  in  a  way  that  increases 
the  interface  area  between  the  fluids  and,  therefore,  allows  enhanced  mixing.  As  the  inflowing  fluid 
reaches  the  first  herringbone,  fluid  at  the  bottom  of  the  channel  is  deflected  toward  the  vertical  side 
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walls.  A  small  portion  is  deflected  one  way  (corresponding  to  the  shorter  segment)  and  a  larger 
portion  is  deflected  the  other  way  (corresponding  to  the  longer  segment).  As  the  fluid  travels 
down  each  segment  of  the  herringbone,  it  finally  reaches  the  side  walls  and  then  is  deflected  up 
these  walls,  over  the  top  of  the  microchannel  (towards  the  center),  and  then  down  towards  the 
bottom.  The  result  is  a  flow  pattern  of  two  convective  rolls  forced  by  the  geometry.  This  is  best 
illustrated  by  velocity  vectors,  as  shown  in  Fig.  9  at  two  cross  sections.  The  larger  convective  roll 
is  located  on  the  side  with  the  longer  segment,  while  the  smaller  convective  roll  is  on  the  side  with 
the  shorter  segment.  At  x  =  600  /im,  the  larger  convective  roll  is  on  the  left  side  of  the  channel, 
and  at  x  —  1540  /xm,  the  larger  roll  is  on  the  right. 

Figure  10  shows  the  velocity  components  at  the  longitudinal  centerline,  at  y  =  45  y m,  which 
is  approximately  20  /im  above  the  herringbone  structures.  The  x-velocity  is  similar  to  that  in  a 
pipe  flow,  with  a  boundary  layer  formed  at  the  entrance.  These  boundary  layers  grow,  merge, 
and  form  a  fully  developed  pipe-flow  profile.  Effects  of  the  herringbone  structures  can  be  seen  in 
the  x-velocity  as  it  increases  slightly  above  herringbones  and  decreases  between  these  structures. 
The  ^/-velocity  increases  when  the  fluid  reaches  a  structure,  and  then  decreases  after  it  passes  the 
structure.  In  other  words,  the  ridges  deflect  the  flow  both  out  and  toward  the  side  walls  and  up, 
and  then  the  the  fluid  flows  down  behind  each  ridge  to  fill  the  space  between  ridges.  Similarly,  the 
velocity  component  changes  direction  depending  on  the  position  of  the  fluid  element  in  relation 
to  the  herringbone  structure.  This  indicates  that  this  is  a  highly  viscous  flow,  driven  by  the 
herringbone  structures;  that  is,  the  cross  sectional  flow  stops  when  the  structure  ends.  As  shown 
in  Fig.  10,  the  magnitude  of  the  ^-velocity  is  approximately  a  factor  of  ten  greater  than  that  of 
the  y-  and  ^-velocity  components. 

We  see  below  that  the  velocity  field  becomes  periodic  in  space  as  the  flow  progresses  down  the 
repeating  herringbone  sequences.  The  first  sequence  is  the  one  that  undergoes  the  most  change, 
due  to  the  proximity  of  the  inflow  boundary.  This  cycle  effectively  sets  up  the  flow  patterns. 
Subsequent  cycles  are  essentially  periodic,  as  shown  in  the  cross  section  at  the  end  of  the  first  and 
second  sequences  in  Fig.  11.  If  we  assume  that  the  velocity  quickly  becomes  periodic,  and  we  know 
that  velocity  field  in  the  second  or  third  sequence,  the  velocity  field  at  any  location  in  a  sequence  of 
herringbones  is  straightforward  to  calculate.  Resolution  tests  with  Ax  =  Ay  =  Az  =  5,  6,  7.5  and 
10  ym  (not  shown)  show  that  the  velocity  field  for  the  5  ym  computational  cell  size  is  converged. 

Given  the  smooth,  resolved,  and  periodic  velocity  field,  we  are  now  able  to  examine  the 
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advection  of  passive  scalars  in  this  field  in  more  sequences  than  were  computed  for  the  velocity 
field.  This  periodicity  of  the  velocity  field  reduces  computational  efforts  significantly.  These 
computations  were  done  in  two  ways:  first,  by  releasing  and  tracking  tracer  particles  that  move  in 
a  Lagrangian  manner  in  the  precomputed  velocity  field;  and  second,  by  using  the  fluid  convection 
algorithm  FCT  to  convect  a  passive-scalar  density  on  the  precomputed  velocity  held.  These  results 
can  then  be  compared  and  assessed. 

3.2  Particle  Advection  on  the  Velocity  Field 

The  first  approach  we  use  to  examine  the  effects  of  the  velocity  held  on  scalar  transport  is  to 
consider  the  behavior  of  a  set  of  Langragian  tracer  particles  released  in  the  how.  Tracer  particles 
of  species  B  are  released  in  the  how  on  the  right  side.  The  equations  solved  for  each  particle  j  are: 

dx^ 

—±=Ui,  i  =  1,  2,  3  and  j  =  1,  N  ,  (8) 

where  x\  and  uj  are  components  of  position  and  velocity,  respectively.  The  boundary  conditions 
are  such  that  when  a  particle  strikes  a  wall,  it  will  stick  to  that  wall  and  is  removed  from  the 
calculation.  No  particle  interacts  with  any  other  particle. 

In  the  calculation  described  here,  we  have  released  9000  particles  on  the  right  side  of  the 
inhow  plane  ( N  =  9000).  This  calculations  uses  the  converged  5.0  /xm-resolution  velocity  held. 
At  the  inhow,  twenty-hve  evenly  spaced  particles  (5  x  5)  were  released  in  each  computational  cell. 
Figures  12-15  show  selected  particle  traces  through  the  four  sequences.  Each  hgure  shows  two 
views:  a  three-dimensional  rendering  that  shows  how  the  particle  moves  across  the  ridges,  and  a 
two-dimensional  projection  of  the  particle  motion  on  the  y  —  z  plane.  These  four  hgures  describe 
very  different  types  of  particle  paths.  Figure  12  shows  the  path  of  a  particle  that  starts  out  on  the 
right  side  and  stays  on  the  right.  It  oscillates  a  little  as  it  moves  downstream,  but  it  stays  within 
a  small  region  of  the  y  —  z  plane.  Figure  13  is  a  very  different  case.  The  particle  starts  out  on 
the  right  side  and,  by  the  end  of  the  first  cycle,  is  far  over  on  the  left  side.  Figure  14  shows  the 
track  of  a  particle  that  stays  in  essentially  the  same  location  on  the  y  —  z  plane  until  early  in  the 
second  sequence,  when  it  abruptly  moves  over  to  the  left  side.  The  final  example,  Figure  15,  shows 
a  particle  that  has  become  trapped  at  a  wall. 

Figure  16  shows  planes  at  the  end  of  each  sequence,  with  white  dots  indicating  where  a  B 
particle  has  intersected  that  plane.  These  figures  result  from  a  one-time  particle  release  and  not 
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a  continuous  stream.  Therefore,  the  time  that  each  particle  intersects  a  particular  plane  could  be 
quite  different.  A  striking  feature  of  the  frames  in  Fig.  16  is  that  the  shapes  of  the  B  regions  on 
the  left  look  as  though  they  have  been  lifted  and  removed  from  the  right  side.  In  fact,  as  we  show 
next,  they  have  been  stretched  and  folded  and  almost  pulled  out.  An  important  point  here  is  that 
we  cannot  really  decipher  the  dynamics  or  the  fluid  evolution  as  it  moves  downstream  from  these 
pictures  alone.  What  we  can  see  is  an  overall  outline  of  the  shape  of  the  regions  of  B. 

There  are  several  sources  of  numerical  error  in  these  particle  traces.  The  first  type  is  the 
residual  error  inherent  in  the  computation  of  the  velocity  field,  which  is  small,  and  to  high  precision 
V  •  v  ~  0.  Other  errors  in  the  particle  computations  could  arise  in  the  specific  velocity  field  used 
for  some  of  the  particles.  All  particles  in  a  particular  computational  cell  were  convected  with  the 
velocity  field  computed  for  that  cell  center,  regardless  of  where  in  the  cell  a  particle  was  located. 
This  is  consistent  with  the  velocity-field  calculation.  Another  option  would  have  been  to  use 
an  interpolated  velocity  field  for  each  location  of  each  particle,  which,  in  turn,  could  introduce 
interpolation  errors. 

3.3  Convection  of  a  Continuous  Variable 

Now  consider  a  different  diagnostic  approach.  Above  we  used  discrete  particles  and  tracked  their 
paths  in  the  flow.  Now  we  consider  a  continuous  representation  and  convect  a  passive-scalar  density 
on  the  same  fixed,  precomputed  velocity  field.  Passive  scalar  densities  were  convected  for  the  5 
/xm  resolution  discussed  above,  as  well  as  for  a  2.5  /xm  and  1.25  /xm  resolution.  The  computations 
at  2.5  /xm  and  1.25  /xm  resolution  were  done  using  a  velocity  field  that  was  interpolated  from  the 
converged  velocity  field  on  the  5  /xm  grid.  Figure  17  shows  snapshots  of  passive  scalars  for  two  full 
sequences  of  herringbone  structures,  computed  at  1.25  /xm  resolution.  The  number  on  each  image 
corresponds  to  the  plane  immediately  after  that  structure.  (That  is,  1  corresponds  to  the  plane 
immediately  after  the  first  herringbone,  and  there  are  12  structures  in  a  sequence.)  These  images 
were  obtained  by  solving  equation 

^  +  V  •  V 4>k  =  o,  k  =  A,B  (9) 

where  </>&  is  the  scalar  density  of  species  A  or  F>,  and  v  has  been  previously  computed.  Because  the 
mass  diffusion  term  had  negligible  effect  in  the  computations,  we  did  not  use  it  in  the  computations 
shown  now.  There  is  a  zero-gradient  boundary  condition  for  A  and  B  at  the  walls  and  at  the 
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outflow.  (This  is  physically  different  from  the  boundary  condition  used  on  tracer  particles  at  walls, 
where  a  tracer  particle  that  hits  a  wall  was  stuck  to  that  wall  and  did  not  move  any  further.)  At 
the  inlet,  there  is  a  continuous  inflow  of  B  on  the  right  side.  The  multidimensionality  was  handled 
by  direction  splitting  the  convection  algorithm  (Oran  &  Boris  2001). 

For  an  incompressible  flow,  the  total  flux  in  each  (y,  2:) -plane, 

/  /  (j)kudydz 

J  J  ( y,z )  plane 

is  conserved  in  the  x  direction.  Because  the  value  of  (j)k  in  each  computational  cell  on  the  inflow 
plane  is  either  0  or  1,  the  values  of  cf)k  in  every  other  downstream  cell  in  the  domain  are  limited, 
0  <  4>k  <  1  because  of  the  conservation  of  (j)k  expressed  by  Eq.  (9).  Values  of  (j)k  can  be  less  than 
1  in  a  computational  cell  because  a  fluid  element  could  be  distorted  by  the  flow  so  that  at  a  later 
time,  its  (constant  but  distorted)  volume  covers  portions  of  several  adjacent  computational  cells. 
When  that  happens,  we  say  that  the  cell  is  mixed  and  the  density  of  is  less  than  1.  Checking 
conservation  and  the  bounds  of  (j)k  is  a  standard  way  to  check  the  accuracy  of  this  computation. 

At  the  inflow  (marked  Start)  shown  in  Fig.  17,  red  fluid,  </)#,  enters  on  the  right  side,  and 
purple  fluid,  enters  on  the  left  side.  For  the  first  six  herringbones,  the  larger  convective  roll 
forms  on  the  left  side.  Therefore,  purple  fluid  is  pushed  towards  the  right  near  the  top  of  the 
microchannel,  while  red  fluid  is  pulled  towards  the  left  near  the  bottom.  When  the  herringbone 
pattern  is  reversed  (structures  7  through  12),  the  red  fluid  is  pushed  upwards  and  over  to  the  left 
at  the  top,  while  purple  fluid  is  pulled  towards  the  right  near  the  bottom.  At  the  end  of  the  first 
sequence  of  herringbones  (after  structure  12),  some  of  the  red  fluid  has  traveled  to  the  left  side, 
and  some  of  the  purple  fluid  has  travelled  to  the  right  side.  This  flow  pattern  repeats  itself  during 
the  second  full  cycle  of  structures,  herringbones  13  -  24. 

The  convection  of  passive  scalar  densities  at  5  /i m,  2.5  fi m  and  1.25  /im  resolutions  were  also 
carried  out  for  four  full  sequences.  The  periodic  feature  of  the  velocity  field  was  used  to  obtain  the 
velocity  field  of  the  third  and  fourth  sequence.  Computing  four  full  sequences  at  2.5  /im  resolution 
required  3200  x  36  x  80  computational  cells,  and  at  the  1.25  /xm  resolution  required  6400  x  72  x  160 
computational  cells.  These  computations  were  done  with  a  fully  parallelized  version  of  a  fluid 
convection  code,  and  required  approximately  10  hours  of  wall  clock  time  on  40  processors  of  an 
SGI  Origin  3800  for  the  2.5  /x m  resolution  and  22  hours  of  wall  clock  time  on  160  processors  for 
the  1.25  /xm  resolution. 
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Figure  18  compares  the  computations  at  the  end  of  four  sequences  for  the  5  fi m,  2.5  /xm  and 
1.25  /xm  resolutions.  For  any  fixed  resolution  (looking  down  any  column  in  Fig.  18),  we  see  that  the 
structure  of  4>b  becomes  more  and  more  complex  as  we  move  downstream  and  pass  more  sequences 
of  herringbones.  The  structure  of  </>#  eventually  evolves  into  the  types  of  striations  seen  in  the 
experiment.  This  behavior  has  been  called  “chaotic  advection,”  described  as  the  complex  behavior 
shown  by  a  passive  scalar  (a  fluid  particle  or  element,  or  a  passively  advected  quantity  such  as 
temperature  or  concentration  of  a  second  tracer  fluid)  when  it  is  driven  by  Lagrangian  dynamics 
in  the  flow  (Cartwright  et  a i.  1999). 

Now  consider  the  plane  at  the  end  of  sequence  1  for  the  three  resolutions.  (This  means  looking 
at  the  top  row  of  the  three  columns.)  From  these  figures,  we  can  see  the  shapes  of  regions  that 
evolve  as  the  accuracy  increases.  This  can  also  be  said  of  the  planes  of  at  the  end  of  sequences  2, 
3,  and  4.  Lower  resolutions  are  not  adequate  as  the  striations  become  finer  and  finer,  but  they  do 
show  how  the  structures  will  evolve  at  downstream  locations. 

The  right-most  column  of  Fig.  18  shows  the  end  of  the  first  four  sequences  for  the  most 
resolved  computation  of  </)#,  at  1.25  /im  resolution.  Scalar  B  particles  moved  across  the  channel, 
from  the  right  to  the  left  side,  and  this  material  appears  as  striations  in  the  flow.  What  is 
particularly  interesting  about  this  flow  is  how  misleading  it  can  be  to  look  at  these  planes  alone, 
without  knowing  how  the  flow  has  evolved  as  it  moves  downstream.  To  make  this  point,  we  have 
reproduced  these  cross  sections  in  Fig.  19,  and  labeled  the  regions  of  B  on  the  left  side  as  a,  5,  c, 
and  d.  For  example,  one  interpretation  might  be  that  material  in  a,  shown  at  the  end  of  sequence 
1,  is  the  same  material  that  is  shown  at  5,  at  the  end  of  sequence  2.  We  have  shown  that  this  is 
not  true  by  looking  at  intermediate  pictures  (such  as  in  Fig.  17)  and  an  animation  of  the  evolution 
of  the  flow  as  it  moves  downstream.  These  show  that  material  a  at  the  end  of  each  sequence  is 
essentially  the  same  material  marked  a  at  the  end  of  earlier  sequences.  And  the  same  is  true  of 
5,  c,  and  d.  That  is,  as  new  islands  of  material  form  on  the  left  hand  side,  the  existing  islands  move 
up  and  towards  the  center. 

There  are  several  sources  of  error  in  the  computations  shown  in  Figs.  17  and  18.  First,  there 
is  the  residual  error  in  the  computation  of  the  velocity  field,  but  we  know  that  this  is  negligibly 
small.  Other  errors  arise  in  the  actual  solution  of  Eq.  (9).  There  is  some  error  introduced  in  the 
solution  of  4>b  by  numerical  diffusion,  even  though  FCT  is  a  high-order  algorithm.  In  the  cases  of 
the  computations  on  finer  resolution  grids,  shown  in  Fig.  18,  there  is  also  some  error  introduced 


11 


by  the  interpolation  of  the  velocity  field  onto  finer  grids. 


3.4  Comparison  of  Particle  Tracks  and  Continuous  Convection 

Figures  16  and  19  can  now  be  compared,  and  we  can  use  the  intermediate  results  of  Figs.  12-15 
and  Fig.  17  to  provide  additional  information.  In  addition,  we  have  an  animation  that  gives  us 
details  of  how  the  flow  and  scalars  change  as  we  step  through  the  computational  domain.  (This  is 
available  on  request.) 

There  are  strong  similarities  in  the  structures  in  Figs.  16  and  19,  which  is  both  encouraging 
and  useful.  Both  approaches  show  the  striations  forming  as  the  tracer  particles  or  scalar  variable 
moves  through  various  sequences.  If  either  Fig.  16  or  19  were  all  that  we  needed,  we  could  obtain 
it  in  a  much  less  expensive  way  by  doing  particle  traces,  rather  than  by  convecting  the  scalar  B. 
Differences  are  apparent  around  the  walls,  and  these  differences  are  consistent  with  the  differences 
in  boundary  conditions  in  the  two  types  of  calculations.  That  is,  in  the  particles  traces,  those 
particles  that  approach  a  wall  are  essentially  trapped  there,  whereas  they  can  eventually  escape  in 
the  scalar  convection  calculation. 

A  striking  feature  of  the  frames  in  Figs.  16  and  19  is  that  the  shapes  of  the  B  regions  on  the 
left  look  as  though  they  have  been  lifted  and  removed  from  the  right  side.  The  history  of  these 
islands  of  £>,  however,  and  not  just  their  final  shapes,  are  critically  important  in  evaluating  the 
performance  of  the  mixer,  particularly  when  the  component  is  applied  to  real  mixing  problems 
with  finite  mass  diffusion.  The  question  is  whether  B  in  these  islands  has  always  been  near  the 
A  —  B  interface  or  if  this  B  material  is  close  to  A  only  after  each  island  is  formed.  If  these  islands 
consist  of  material  that  had  previously  been  adjacent  to  the  right  wall,  then  this  material  B  was 
isolated  from  material  A  until  the  mixer  pulled  it  to  the  left  and  placed  it  in  close  proximity  to  A. 
This  displacement  of  fluid  B  would  produce  huge  concentration  gradients  in  the  diffusing  species 
and  rapid  mass  transfer  between  the  fluids  at  the  islands.  If,  in  contrast,  the  material  in  the  islands 
originates  near  the  A  —  B  interface,  then  mass  diffusion  would  have  already  smoothed  the  gradient 
between  the  two  fluids,  and  the  mass  transfer  at  the  island  would  be  reduced. 

To  understand  how  these  regions  form,  we  examined  the  dynamics  of  this  flow.  This  was  done 
by  looking  at  the  evolution  of  B  as  the  flow  moved  over  the  sequences  of  structures.  We  found  that 
the  material  marked  a  in  Fig.  19  was  originally  on  the  interface  between  A  and  B  at  the  inflow. 
Then  as  the  material  passed  over  later  sequences,  regions  5,  c,  and  d  were,  in  turn,  moved  from  the 
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A  —  B  interface.  These  regions  of  material  were  stretched,  folded,  and  ultimately  pinched  from  the 
larger  B  region  along  the  A  —  B  interface.  This  evolution  can  also  be  seen  in  Fig.  17.  Starting 
at  herringbone  1,  red  fluid  from  the  A  —  B  interface  is  moved  over  to  the  left  along  the  bottom 
of  the  channel  for  the  first  six  herringbones.  Then,  when  purple  fluid  is  moved  over  to  the  right 
along  the  bottom  of  the  channel  for  the  second  six  herringbones,  the  red  fluid  on  the  left  (from 
the  A  —  B  interface)  is  pinched  off  from  the  larger  B  region,  resulting  in  the  island  of  red  fluid  on 
the  left  hand  side  at  the  end  of  the  first  sequence  (Herringbone  12).  This  movement  of  material 
is  confirmed  by  selected  particle  tracks,  but  the  advection  of  the  continuous  scalar  concentration 
demonstrates  this  evolution  much  more  clearly. 


4.  Comments  and  Conclusions 

Mixing,  due  to  enhanced  stretching,  bending,  and  folding  of  material  surfaces  in  a  flow,  has  been 
studied  extensively  in  systems  that  show  chaotic  advection.  Mixing,  however,  is  only  one  part  of 
the  problem.  We  are  also  interested  generally  in  flow  control,  and  that  includes  mixing,  separation, 
delivery,  and  transport.  The  development  of  concepts  of  nonlinear  dynamics  can  be  helpful  in  all  of 
these  areas.  In  terms  of  nonlinear  dynamics,  a  fixed  point  in  the  system  is  one  in  which  a  particle 
starts  out  at  some  point  (x0,y0,z0)  and  remains  for  all  time  t  at  (x,y0,z0).  The  particle  travels  in 
a  straight  line  through  the  system.  A  periodic  point  is  one  in  which  a  particle  released  will  return 
to  its  initial  position,  again  and  again,  after  a  fixed  time  or  space  interval.  At  intermediate  times, 
it  is  elsewhere,  and  could  be  anywhere.  What  is  important  is  that  you  know  when  to  find  it  at  the 
specified  point,  and  in  our  case,  we  know  where  it  has  been.  Sometimes,  there  are  periodic  points 
but  their  period  is  longer  —  it  takes  more  cycles  for  the  particle  to  return  to  its  initial  point. 

Passive  scalars,  as  those  shown  in  Figs.  12-19,  are  particles  that  are  not  large  enough  to 
perturb  the  flow,  but  large  enough  not  to  diffuse  (that  is,  there  is  no  Brownian  motion).  As 
such,  they  only  approximate  the  motions  of  real,  and  often  relatively  large  particles  of  interest  in 
a  background  flow.  Our  velocity  field  is  steady,  but  periodic  in  space.  Therefore,  we  could,  in 
principle,  look  for  fixed  points  and  periodic  points  in  the  flow.  Figure  16  is  a  Poincare  map,  a 
reduction  in  the  dimension  of  the  system  found  by  using  a  section  of  the  flow  (in  time  or  space) 
to  produce  a  map.  The  dynamic  behavior  of  this  map  is  sometimes  easier  to  study  than  the  full 
system,  and  sometimes  gives  some  useful  insights  into  the  system.  For  example,  the  particle  shown 
in  Fig.  12  is  close  to  a  fixed  point.  Other  points  might  be  close  to  periodic  points.  In  general, 
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however,  we  seem  to  learn  more  about  the  dynamics  of  the  system  and  the  physical  mechanisms 
of  the  flow  by  focusing  on  the  fluid  elements  and  particle  motions  in  three  dimensions. 

In  the  experiment  (Fig.  2),  a  laser  aimed  at  the  liquid  containing  the  light-sensitive  dye 
provided  visualization  of  the  mixing  between  the  two  streams.  Note  that  the  visualization  was 
not  calibrated  -  we  do  not  know  exactly  how  much  mixing  corresponded  to  a  certain  intensity 
of  light.  The  scale  could  have  been  linear,  or  simply  a  sharp  cut  off,  or  some  highly  nonlinear 
function.  In  addition,  the  laser  experiments  were  affected  by  light  reflection  from  the  walls.  This 
would  make  the  signal  appear  stronger  near  the  walls.  This  was  indicated  to  some  extent  in  the 
experimental  paper,  which  used  only  an  inner  portion  of  the  image  (the  part  enclosed  by  the  dashed 
lines  in  Fig.  2)  for  the  analysis.  Finally,  at  a  relatively  late  stage  in  our  work,  we  learned  that 
the  experimental  herringbones  are  actually  grooves  in  the  channel,  not  ridges  protruding  into  the 
channel  flow  (Stroock  et  a 1.  (2002);  Stroock  &  McGraw  (2004)),  as  was  simulated  here.  However, 
given  this,  and  comparing  Figs.  2  and  18,  we  can  conclude  that  grooves  and  ridges  produce  similar 
qualitative  effects  at  the  scales  in  which  we  are  looking. 

Numerical  resolution  is  an  important  issue  here.  We  really  cannot  resolve  the  striations  that 
form  after  multiple  sequences  in  any  reasonable  amount  of  computer  time.  Then  we  ask:  do  we 
really  need  to  do  this?  It  is  unlikely  that  practical  devices  will  need  more  than  a  few  sequences,  and 
recent  studies  (Howell  et  a 1.  2005)  have  shown  that  by  varying  the  pattern  ridges  or  grooves,  we 
can  achieve  mixing  and  flow  control  in  much  shorter  distances  than  15  cycles  or  even  four  cycles. 
In  any  case,  we  can  compute  v  to  high  enough  accuracy,  and  this  is  really  the  beginning  of  the 
problem.  We  can  use  this  v  to  study  the  dynamics,  as  long  as  the  particle-fluid  interactions  are 
small  and  do  not  affect  v.  When  these  interactions  do  affect  v,  the  computations  of  v  and  the 
particle  behavior  must  be  done  simultaneously  and  interactively,  more  as  it  is  done  in  multispecies 
or  multiphase  fluid  dynamcis. 

A  major  conclusion,  then,  is  that  we  were  really  not  able  to  analyze  and  fully  understand  the 
flow  from  the  particle  traces  alone.  Without  considering  the  details  of  the  flow  evolution,  which 
were  obtained  from  the  passive  scalar  advection,  we  could  not  understand  how  certain  parts  of 
scalar  B  were  moved  across  the  channel.  These  clarified  what  we  saw  in  the  particle  traces.  Both 
approaches,  however,  provided  valuable  information  and  were  used  to  back  up  or  substantiate  the 
conclusions. 
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Figure  1.  Schematic  of  a  microchannel  with  herringbone  ridges.  The  first  two  sequences  are  shown. 
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Sequence 


Figure  2.  Confocal  micrograph  images  taken  at  the  end  of  the  first  five  sequences  and  the  last  sequence 
(> Stroock  et  a 1.  (2002)). 
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Figure  3.  Computed  x-velocity,  z-velocity,  pressure,  and  concentration  of  B  in  the  x-z  plane  at  y  =  45  pm 
for  the  unobstructed  microchannel. 
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Figure  4.  Concentrations  of  scalars  A  and  B  at  the  inflow  (x  =  0  pm)  and  exit  (4000  |Llm)  planes, 
respectively,  for  the  unobstructed  microchannel. 
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Figure  6.  Example  of  the  tests  performed  to  evaluate  the  effects  of  staircasing  in  the  computational 
grid.  A  triangle-based  obstacle  is  placed  in  the  flow,  with  staircasing  in  the  z-direction. 
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Figure  7.  Geometry,  computed  velocities,  and  pressure  for  the  obstacle  shown  in  Figure  6.  Quantities 
are  shown  in  the  x  -  z  plane  at  y  =  25  pm  for  the  two  staircasing  resolutions. 
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Figure  8.  Profiles  along  the  x-axis,  at  y  =  25  |0,m,  z  =  100  |0,m  for  the  two  staircasing  resolutions.  The  flat 
region  in  each  profile,  where  215  <  x<  250  |0,m,  corresponds  to  the  obstacle  location. 


25 


Figure  9.  Velocity  vectors  in  the  y  -  z  plane  showing  the  flow  pattern  of  geometry-induced  convective 
rolls  in  the  cross  section  of  the  microchannel  shown  in  Fig.  1  at  two  x-locations,  600  |tm  and  1540  pm. 


T 

z 


0.21  —  2.28  Vx  (crn/s) 


-0.34  _  0.34  Vz  (cm/s) 


Figure  10.  Velocity  components  at  the  longitudinal  centerline  y  =  45  |0,m  for  the  microchannel  shown  in  Fig  1. 


End  of  Sequence  1 


End  of  Sequence  2 


Figure  11.  Velocity  field  at  the  end  of  sequences  1  and  2,  showing  that  it  becomes  essentially  periodic 
at  the  end  of  each  sequence. 
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Figure  12.  Track  of  a  particle  that  maintains  essentially  the  same  y  -  z  position  (i.e.,  is  essentially 
unaffected  by  the  ridges)  as  it  moves  down  the  channel,  (a)  Projection  on  the  particle  path 
in  the  y  -  z  plane,  (b)  Three-dimensional  view  of  the  particle  motion. 
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Figure  13.  Track  of  a  particle  that  is  shifted  during  the  first  sequence  to  the  other  side  as  it 
moves  down  the  channel,  (a)  Projection  on  the  particle  path  in  the  y  -  z  plane.  The  numbers 
in  the  particle  projection  correspond  to  the  end  of  a  sequence,  (b)  Three-dimensional  view. 
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Figure  14.  Track  of  a  particle  that  is  shifted  during  the  second  sequence  to  the  other  side  as  it  moves  down  the 
channel,  (a)  Projection  of  the  particle  path  in  the  y  -  z  plane,  (b)  Three-dimensional  view. 
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Figure  15.  Track  of  a  particle  that  is  essentially  stopped  at  a  wall,  (a)  Projection  on  the  particle  path 
in  the  y  -  z  plane,  (b)  Three-dimensional  view. 
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Figure  16.  Poincare  map  showing  the  locations  where  particles  B  crossed  the  y  -  z  plane  after  each  sequence. 
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Figure  17.  Concentration  of  passive  scalars  A  (purple)  and  B  (red),  shown  in  the  y  -  z  plane,  after  each  herringbone  in 
the  first  two  sequences.  Computation  at  1.25  pm  resolution.  Numbers  in  upper  right  corner  of  each  frame  (from  1  to  24) 
indicate  locations  as  shown  on  the  schematic  on  the  left  side  of  the  figure. 
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Figure  18.  Concentration  of  passive  scalars  A  (purple)  and  B  (red)  on  the  y  -  z  plane  at  the  end  of  each  sequence,  computed 
at  three  resolutions,  using  the  velocity  computed  at  5  jlm  resolution. 
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Figure  19.  Concentration  of  passive  scalars  A  (purple)  and  B  (red)  on  the  y  -  z  plane  at  the  end  of  each  sequence,  computed 
at  the  highest  resolution  (1.25  |0,m).  Folding  and  stretching  of  fluid  elements  results  in  B  being  moved  from  the  right  hand 
side  to  the  left  hand  side  of  the  mixer.  This  results  in  the  formation  of  island  a  on  the  left  hand  side  at  the  end  of  first 
sequence.  At  the  end  of  the  second  sequence,  island  b  has  formed  on  the  left  hand  side,  while  island  a  has  moved  upwards 
and  towards  the  center.  Similar  sequences  of  bending  of  folding  produce  islands  c  and  d. 
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